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Abstract 

The problem of the focal region fields scattered by an arbitrary surface reflector under 
uniform plane wave illumination is solved. The Physical Optics (PO) approximation is used 
to calculate the current induced on the reflector. The surface of the reflector is described 
by a number of triangular domain-wise 5th degree bivariate polynomials. A 2-dimensional 
Gaussian quadrature is employed to numerically evaluate the integral expressions of the 
scattered fields. No Freshnel or Fraunhofer zone approximations are made. The relation 
of the focal fields problem to surface compensation techniques and other applications are 
mentioned. Several examples of distorted parabolic reflectors are presented. The computer 
code developed is included, together with instructions on its usage. 
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Introduction 


Far field electromagnetic scattering from conducting surfaces of arbitrary shape has 
been studied by several researchers [ 1—6] . To our knowledge, however, little or no attention 
has been given to the near field scattering from arbitrary surfaces. The analysis of this 
problem may be useful to the efforts of compensation of reflector surface distortions. To 
illustrate this, let us consider the following situation. A plane wave, E , is incident on 
an arbitrary reflector; the scattered fields, E * , in a region, n, in the neighborhood of the 
reflector are calculated by the method outlined in the present paper. Using the reciprocity 
theorem yields, if in fl the fields E 3 are considered to be incident on the reflector, then the 
fields scattered by it will be E* . For antenna performance considerations E l is required 
to meet certain specifications. These specifications can be easily employed in the context 
of the method presented here and yield the E 3 fields in 0. If 0 is chosen to be the 
region occupied by the feed, then E 3 is related to the excitation. Thus a scheme can be 
developed to provide the necessary excitation that renders desired far fields, E\ from a 
reflector antenna of arbitrary shape. In addition, focal field analysis may be applied to the 
problem of scattering by objects in the neighborhood of the feed and/or on the feed itself. 

In the present paper, we develop a straightforward method to solve for the scattered 
fields (with emphasis in the near field zone) when a plane wave is incident on a reflector 
with arbitrary surface. The currents induced on the reflector are calculated by the Physical 
Optics (PO) approximation (K = 2n x H *). Subsequently, these currents are used as 
sources of the scattered fields and explicit integral expressions of the latter in terms of K 
are shown. These integrals expressions involve no approximation due to the position of 
the observation points with respect to the sources and, therefore, are valid for all zones 
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(near, Fresnel and Fraunhofer). Since we deal with rapidly oscillating integrands, accurate 
integration techniques are required and employed in the present work. 

The novelty of our method consists of: 

1. Triangular discretization of the reflector aperture by a max-min criterion based 
algorithm. 

2. Reflector surface representation via interpolation of several fifth degree bivariate 
polynomials. 

3. Development of a systematic, 2-dimensional, Gaussian quadrature calculation of the 
integrals involved in the scattered field expressions. 

In this report, we examine the effects of various surface distortions on the focal region 
fields of parabolic reflectors. The performance of distorted reflectors with frequency is also 
discussed. 


Theory 

Let us consider a plane electromagnetic wave incident on an arbitrarily shaped and 
perfectly conducting surface, E, as in Figure 1. Let us also consider a Cartesian coordinate 
system Oxyz such that the magnetic field of the incident wave is along the y axis and the 
wave vector is along the negative z axis. The surface of the reflector in this coordinate 
system is described by the function 

z = g(x,y). (1) 

The fields of the incident wave are given by 

H = yH 0 e ]/3z (2) 
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and 


E = —xrjHoe 3 


where 




and (3 is the wavenumber. x,y and z are the unit vectors along the axis of the 0 xyz 
coordinate system and 17 is the intrinsic impedance of the medium. According to the 
physical optics ( PO ) approximation, a surface current density K is induced on the reflector. 
This current density is the source of the scattered fields and is equal to 


K{f) = 2 n(f) x H{f) ; FcE 


where n(r) is the unit vector normal to the surface of the reflector. Using (1) we can obtain 


the expression 


_ g- Vgfc.y) 

\z - V<7(x, y)| 


Elaborating on Equation (6) and using Equation (2) and (5), the surface current density 
K(r) is derived to be 


K(r) = 2 H 0 dx - e 3 


The vector potential due to this current is 


— If— e -jP\ r - f '\ 

A (f) = ~ / K(r') - — - d 2 r’ 


The magnetic and electric scattered fields are related to A(r) by 


H{f) = V x A{f) 
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and 


Elf) = ^-V x V x Air). (10) 

juje 

Expanding Equation (10) we obtain the following integral expression of the electric field 
in terms of the current density. 


E(f) 


= -±-f 

4irjue 


e-W R 


47T ju>e Jz R 

1 - /3 2 R 2 + j(3R 

R 2 


K(f') • (f - f') 


K{f') dV. 


3 + 3 j(3R-f3 2 R' 

R 4 


(f - f') 


( 11 ) 


where e is the dielectric permittivity of the environment and R = f — f'. Equation (11) is 
in agreement with an alternative expression given by Silver [7]. The integration is carried 
over the illuminated portion of the reflector and d 2 r' is the area element on E. If o is the 
projection of E on the x — y plane, the integral of Equation (11) can be performed on o due 
to the relation dictated by Equation (1). Differential geometry considerations [2,8] yield 
the following relation between the area elements on E and o,d 2 r' and dx'dy' respectively. 


d 2 r' = 


1 + 


dg{z\y') 

dx' 


+ 


dg(x',y') 

dy' 


dx'dy’ . 


( 12 ) 


Combining Equation (7) and (12) we obtain 


E(f , )d 2 r' = 2H 0 e jl3 *‘ 



~ dg{x',y') ' 

dx' 


(13) 


Substituting Equation (13) into Equation (11) we get 

e -jf3(R~z’) 


E(r) = - / 

47r;u>6 J a 


R 


-(x - x’) - (z - z')^, 


2 d2 


3 + 3 jdR -P 2 R 

R 4 


.(f-f’) 


\ 1 -0 2 R 2 +j0R / A * dg 


+ 


R 2 


x + 2 — )dx dy . 

ox 


(14) 
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The expressions of the scattered electric field in Equation (14) has been obtained by 
making only one approximation, namely the PO approximation (cf. Equation (5)). Since 
we are primarily interested in the focal region fields, the usual and convenient far field 
approximations cannot be used here. Neither can the Fresnel region approximations since 
we would like to treat geometries for which the focal point is in the near zone. Since we 
want to analyze arbitrary surfaces, a good way to represent the reflector is by a number 
of “target points.” In an experimental version of this problem, target points would be 
small dots of an appropriate material placed on the surface of the reflector. We assume 
that the coordinates of these points can be measured by some technique (photogrammetic 
measurements). The projection of these points on the x — y plane provides a natural 
discretization of the domain of integration a (see Figure 2). We employ a max min criterion 
algorithm [9,10] to define triangles with vertices which are the projections of the target 
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points on a. Since the collection of all triangles constitutes a polygonal approximation 
to cr, the integrals over a can be found by simply calculating the contribution from each 
triangular domain and subsequently adding them up. For an analytical representation of 
the surface of the reflector,, the same program mentioned above is used to interpolate a fifth 
degree bivariate polynomial in each triangular domain. This polynomial representation is 
used to estimate the values of z' and -§^7 on the surface of the reflector. Needless to say, the 
larger the number of target points and, thus, triangles on a, the better the representation of 
E. The integral of Equation (14) is not only vectorial but, it is complex as well. Therefore, 
numerically we need to perform six integrations in order to find the real and imaginary 
parts of all three components of the electric field at an arbitrary point f = xx + yy + zz. 
The functions to be integrated can be very oscillatory in nature due to the phase term 
exp(— j/3(r — z')) and also due to the presence of the dg/dx' term. As is well known [11,12] 
there are no satisfactory quadrature formulae to integrate functions on triangular domains. 
However, a series of mappings can be employed to transform a triangle into a square (see 
Figure 3). There, a simple extension of the one-dimensional Gaussian Quadrature is 
available. A brief outline of this method follows. 

By means of a linear transformation, every triangle of the x-y plane can be mapped 
onto a specific triangle in another plane, u-v. Consider an arbitrary triangle in the x — y 
plane with vertices (x^yj), ( 12 , 3 / 2 ) an ^ ( x 3> 2 / 3 )* This arbitrary triangle can be mapped 
onto a basic triangle in another coordinate system, u — v, by the transformation: 


x = xi + au + bv 

(18a) 

y = yi + cu + dv 

(186) 


where 

a = x 2 — xi (19) 
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6 = x 3 - x\ 


( 20 ) 

( 21 ) 


c = y2 - y\ 


and 


d = yz - y\ 


( 22 ) 


The Jacobian of mapping (18) is 


( ZiIL) = 

a 6 

\u,vj \ 

c d 


= [ad - be) = (x 2 - xi)(y 3 - y i ) - (*3 ~ *i)(ya ~ Vi). 


(23) 


The inverse transformation is 


u = [ d{x - ii) - b(y - yi)]/(ad - be) (24a) 

v = [ c(x - ii) + a(y - yi )]/ (ad - 6c). (246) 


According to (18), the vertices A = ( xi,yi)B — (x 2 ,y 2 ) and C = (13,2/3) are mapped 
onto the points [0,0], [1,0] and [0,1], respectively, in the u — v coordinate system. These 
points define what we call the basic triangle in the u — v system. As mentioned above, 
there are no satisfactory quadrature formulae for integration on triangular domains for 


rapidly oscillatory integrands. The basic triangle in the u — v system can be mapped onto 
a square domain by means of the transformation 

1 + r 


u — 


The Jacobian of mapping (25) is 


J 


v=[ 

'1 -r\ 

4 J 

(1 + . 

5) is 



( _ 

1/2 

0 

V r ,s) 

14-3 

4 

1 — r 
4 


(25a) 

(256) 


1 - r 


(26) 
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According to (25), the sides v = 0 and v = 1 — u of the basic triangle are mapped onto the 
lines 5 = — 1 and 5 = +1 respectively, in the r — s coordinate system. Similarly, the side 
u = 0 is mapped onto the line r — — 1 and, the point B is mapped onto the line r — 1. In 
Figure 2, we show the series of mappings described above. 

An integral over D xy can be numerically calculated using a Gaussian quadrature 
formula on the square domain D ra provided we represent the integrand appropriately. 
Thus, if f(x,y) is a function defined over the triangular domain D xy and I is the integral 

1= J J f{x,y)dxdy, (27) 

D Z y 

we know that 

/ = J J f(x(u,v),y(u,v))J^^ S jdudv (28) 

and also 

I = j J /^x^u(r,s)^,y^u(r,s),u(r,s)^J J (7^) dr ds. (29) 


For our purposes, Equation (29) essentially looks like 


/= J J T(r, s)drd$. 


(30) 


where 


r(r, S )=/(*,v)j(^)/(^) (31) 

An obvious 2-dimensional extension of the well-established 1-dimensional quadrature 


formulae is used here. That is 


and, finally 


N N 

' = EEwh.<ii' (32) 

k=i j=i 

In our analysis r k and Sj ( k,j = 1,2,..., TV) are the roots of the N </l order Legendre 
polynomial, P jv, and W k {k — 1,2,..., TV) are the appropriate weights for Gaussian 
quadrature, given by 

' (1 -r 2 F )P' N (rlY 

From (32) it is clear that the values of T are required at the points {r k ,Sj) of the r — s 
coordinate system. Using the transformations (25) and (18), the corresponding points 
(x k j,y k j) of the x — y coordinate system can be found. It is clear that 



f{ x kj jt/kj) — F{ r ki s j)' (34) 

There are TV 2 points in each triangle. The value of TV required to accurately calculate the 
integrals depends on the frequency and the surface of the reflector. It is found that if the 
reflector has extensive distortions, a large TV is required. 

One subtle point needs to be mentioned here. Since the square domain D ra and the 
basic triangle domain D uv are fixed and kept the same for all x - y triangles, the limits of 
integration in (28) and (29) are desired to stay fixed. As a consequence, A, B, and C must 
be in a counter-clockwise sense. It is important that the vertices of all the D xy triangles be 
arranged in the same rotational sense. If instead of the counter-clockwise sense presented 
here, a clockwise sense were chosen, then the Jacobian J would need a sign reversal 

to ensure the correct result. 

To recapitulate, the projection, <r, of the reflector on the x — y plane is divided in a 
number of triangles. The integrals of Equations (15)— (17) are decomposed into summations 
of integrals over the aforementioned triangles. For each triangle, the points (x kj ,y kj are 
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found by applying (25) and (18). Then, the quadrature formula (32) is used to evaluate 
the contribution of each triangle to the total integral. Finally, all these contributions are 
added together to yield approximate values for E sx ,E sy and E sz . 

Oblique Incidence 

The generic geometry presented in the second section can be used to solve the problem 
of the scattered fields of a given reflector under oblique incidence. To demonstrate this 
assertion, consider a coordinate system Oxyz, attached to the reflector. The relative 
orientation of the reflector system, 0 xyz, and the incident wave system, Oxyz, is described 
by the Euler angles 9,<f> and i/> as shown in Figure 4. 

To solve for the scattered fields under these circumstances, we make use of the fact 
that g(x,y) in Equation (1) is arbitrary. If the reflector in its own coordinate system is 
described by z = y(x,y), then we can express this reflector in Oxyz and then use the 
analysis of Section 2. To do this we perform three rotations. First a positive rotation 
about Oz by <f>. This rotation moves 0 y to a new position 0 y . The second rotation is 
about Oy by +6. This rotation aligns the optical axis of the reflector with Oz. Finally, 
a third rotation about Oz by an amount ip is required to align Oy with Oy and thus take 
care of the polarization. 

The sequence of rotations is shown in Figure 5. To demonstrate the use of these 
rotations, consider the following special cases. 

1. Incident wave E on the Oxy plane and H field perpendicular to the 0 xz plane (see 

Figure 6). Thus, <p = 0 and ip = 0. 
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2. Propagation along the negative 0 2 axis (see Figure 7). Thus, 8 — 0, 0 = 0, and 
0 = 0. Please note that you need to rotate the 0 xy z system to make it align 
with the Oxyz system. Furthermore, the values of 0,0, and 0 have to be positive for 
counterclockwise rotations and negative for clockwise ones. 

Numerical Examples 

In an experimental application of the analysis presented above, the coordinates of the 
target points would be measured and provided as data to the developed computer code. 
Since this report does not involve experimental studies, an analytical method is needed to 
generate the coordinates of the target points. First we distribute the projections of the 
target points onto the Oxy plane. We do this by locating all of these points on circles of 
successively increasing radius up to a given value. The number of points on each circle 
also increases as we move outwards. This scheme of discretization is suitable for circular 
aperture reflectors and is shown in Figure 8. In this discretization example, we have used 
four circles with radius and points as follows. 


Circle # 

Radius 

Points on Circle 

0 

0 

1 

1 

1 

np = 3 

2 

2 

2np = 6 

3 

j 

3 

2(2np) = 12 


This scheme provides us with the x and y coordinates of the target points. To complete the 
description of the surface of the reflector, a function is needed to generate the 2 coordinates 
of the targets. For all the examples presented in this section, the target points lay on a 
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surface generated by the function 


2 


rv 2 2 

x j- y 

4F 


+ A\t 


A^Kx— As)' + (y — A^)') 


~2 D 2 
+ y < T 


(35) 


Clearly, Equation (35) represents an ideal paraboloid of focal length, F, that is perturbed 
by a Gaussian “bump” of height, Aj, steepness, A 2 , located at the point (x,y) = (A3, A4). 
For Ai = 0 Equation (35) represents an ideal parabolic reflector of diameter, D. We are 
primarily interested in the fields of the focal region. The observation points in this section 
will be assumed to be distributed on a line segment of electrical length /?£ that has its 
origin at the point Of, on the z-axis and is parallel the the 0 xy plane (see Figure 9). 

In Figure 10, the magnitude of the electric field components is shown for an ideal 
paraboloid with diameter, D = 2m, focal length, F = 1.19m, at frequency, / =500 MHz. 
The observation points are defined by <£„(, = 0 and /?£ = 20, i.e., they lay on the 0j,x axis. 
In this example <$> — 9 — tp = 0°, i.e., the coordinate systems Oxyz and 0 xy z are identical. 
That is the incident field is propagating along the negative Oz axis and has its H(E) field 
along the 0y(0x) axis. It is clear from Figure 10 that the cross polar field, E y , is very weak, 
negligible in comparison to E x and E z . In Figure 11, the magnitude of the electric field 
components of a perturbed version of the paraboloid of the previous example is shown. 
The characteristics of the perturbation are A\ =0.1m, A 2 = -3, and ^3 = A4 =0.2m. A 
comparison between Figure 10 and 11 reveals that the surface perturbation has quantitative 
effects on the focal region fields. First, the co-polar component has a maximum decreased 
by 10% from the corresponding value of the ideal surface case. Second, the cross-polar 
level is somewhat increased and third, the E z component is no longer zero at the point 0 

In Figure 12, the magnitude of the normalized electric field components of the same 
distorted reflector used in Figure 11 is shown. The only difference between Figure 11 and 
12 is the frequency. In Figure 11, the frequency is 5 GHz. The observation points are 
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again distributed on the 0{>x axis and /?£ is the same as earlier (/?£ = 20— > t « 20 cm). 
As is expected, because of the higher frequency, the bump on the surface of the reflector 
has more severe effects now than in Figure 11. To illustrate this, we compare the focal 
fields of the distorted reflector with the fields of the ideal reflector at / = 5 GHz. This 
comparison is depicted in Figures 13a-13f. From Figure 12, we see that the cross-polarized 
field component ( E y ) is eight times smaller than the coplanar component or, down by 18 
dB. We also see that, as in Figure 11, the axial field at the point 0*, is not zero when the 
surface has been distorted. Furthermore, we note that the phase of the cross-polar field 
component of the ideal surface parabolic reflector is not significant since the magnitude of 
the cross-polar field is zero (Figure 13b). Because of this fact, the phase of the cross-polar 
component of the ideal reflector is shown to exhibit a rather large numerical inaccuracy 
(Figure 13e). We may also point out that there is a significant increase (from 3.5 to 37) 
of the normalized coplanar component when we go from / = 500 MHz to / = 5 GHz (see 
Figures 10 and 13a). This is in accord to the well known limit of infinite focal fields at the 
Geometric Optics (GO) approximation. 

Exactly the same number of target points has been used for all of the examples 
presented above. Furthermore, the coordinates of these points in the Oiy plane were 
also exactly the same. This enabled us to compare the results of the various test cases 
used. It is important to realize that, in general, the domain-wise 5th degree polynomial 
interpolation only approximates the surface of the reflector. As a consequence, if one uses 
two different sets of x-y coordinates of target points to represent a given reflector surface 
(for example given by Equation (35)), different scattered fields will be obtained under 
the same conditions of illumination. This problem, however, can be solved very easily by 
increasing the number of target points to achieve good polynomial representation of the 
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reflector surface. Numerical verification of this solution has been done, and we have found 
that the higher the frequency, the larger the number of target points required to represent 
the reflector in order to achieve convergence of the scattered fields. 

At / = 35 GHz, and when the propagation is along the optical axis, the focal region 
fields of ideal parabolic receiving antennas obtained by the method of the present report 
have been found in complete agreement with the numerical results obtained by a different 
method (expansion in spherical wave functions) as in [13]. For arbitrary illumination 
and/or reflectors, however, we have not found studies by other researchers to cross-examine 
our numerical results. 

In these cases, the scattered fields at a few observation points are calculated using 
various densities of points in the Gaussian quadrature (this is equivalent to varying N in 
Equation (32)) until convergence is achieved. For a given frequency, we have found that 
the smoother the surface of the reflector (the smaller A\ in Equation (35)), the smaller N 
is required for convergence of the fields. To demonstrate this process, the ideal parabolic 
reflector used earlier ( D = 2m, F = 1.19m) is subject to obliquely incident illumination 
(6 = 5°,4> = 0, = 0). In contrast to the previous examples, now the antenna (0 xyz) 
and the ray (0 xyz) coordinate systems are not identical. The normalized incident field, 
E* / Eq in the two systems is 


E^_ 

E o 


= — xe 


_ _ 


( : S’ 


— z sin 5 



sin 5 ° + z cos 5°) 


(36) 


The observation points are distributed on the negative x-axis, i.e., pob = 180° and Pi = 
20— ► £ ~ 20 cm. The geometry is shown in Figure 14. The magnitude of the normalized 
electric field components in the antenna coordinate system is plotted in Figures 15, 16, 
and 17 for various values of N = ng. From these figures, it is evident that convergence is 


15 



achieved very quickly. Note that the y-component is the smallest of all; about four orders 
of magnitude smaller than the x and z components. Convergence of the y-component is, 
however, still very rapid. In Figure 15, we observe that the scattered field has a maximum 
at a distance xo = —11.2 cm. If we interpret this in degrees from the optical axis, we find 
that the maximum of the scattered field occurs at an angle, 9o, where 

~ 5.38° . (37) 

This was expected, since the angle of incidence, 9, is 9 = 5°. 

As an additional example, consider the previous case when the reflector is distorted 
by a bump located at the origin. The surface of the reflector is described by Equation 
(34) with parameters D — 2m, F — 1.19m, A\ — 0.1m, A 2 = -3, A 3 — 0, and A 4 — 0. 
The frequency is again / = 5 GHz, and the incident illumination is described by Equation 
(36); i.e., 9 = 5 °,<j> = 0,xp = 0. The magnitude of the normalized components of the 
scattered field on the negative x-axis (0o6 = 180°, /?£ = 20) is shown in Figures 18, 19 
and 20 for various N = ng. The convergence is achieved rapidly, especially for the major 
components of the field (E~ and E~). Comparing the undistorted (Figure 15, 16, and 17) 
with the distorted (Figure 18, 19, and 20) we observe that the largest field component, 
suffers a severe reduction (by a factor of 5) from its level at the absence of distortion. 
The E~ component suffers a comparable reduction. Moreover, the level of the smallest 
field component, E~, is increased (by an order of magnitude) due to the distortion. 

The parameter A\ — 0.1m used in the examples of distorted reflector examples is 
a rather extensive distortion, since it represents 10% of the radius of the reflector. In 
practical applications less severe distortions are expected. The higher the frequency, /, 
the larger the angle of incidence, 9, the larger N is required to achieve convergence. 
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Conclusions 


A straightforward method to calculate the magnitude and phase of all three com- 
ponents of the focal region field scattered by an arbitrary reflector under plane wave 
illumination is presented. Discretization of the surface of the reflector and subsequent 
interpolation by fifth degree bivariate polynomials is made. The current induced on the 
reflector is evaluated by the PO approximation. The scattered field is given in an integral 
expression that involves no Fresnel or Fraunhofer zone approximations. A novel approach 
to the problem of numerical integration of the integrals is discussed. This method is based 
on a 2-dimensional Gaussian quadrature and provides a high degree of flexibility since 
the density of the sample points can be varied interactively as to adapt to the accuracy 
needs of the specific application. Several examples of distorted and undistorted parabolic 
reflectors are presented. The behavior of the scattered field with frequency and angle of 
incidence is discussed. 
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Figure 1. Arbitrary Reflector Geometry 
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Figure 2. 


Reflector Discretization 
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Figure 3. Successive mapping of an arbitrary triangle 
onto a basic one and, then, a square 
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Figure 4. Relative orientation of antenna and ray 
coordinate systems . 
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Figure 5 . 


Euler 1 s angles 0, <)> and ip define the 
relative orientation between the antenna 
and the ray coordinate system . 
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Figure 8. Discretization Scheme 
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Figure 10. Normalized E. Observation points on the Ox axis. 

/ = 500 MHz, D = 2m, F = 1.19m, /?£ = 20, 
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Figure 13a. Normalized \E X \ of ideal and distorted reflectors. 
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Figure 14. Oblique Incidence at / = 5 GHz. 0 = 5°, <f> = 0, xJj = 0. 
D = 2m, 01 = 20, 00 6 = F = 1.19m. 
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Figure 15. 


Normalized \Ex\. Oblique incidence (0 = 5° , <p = 0, tp — 0) 
D — 2m, F = 1.19m, A\ = Ai — A 3 — A\ = 0. 
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Figure 16. Normalized \Ey\. Oblique incidence (9 = 5°, <f> = 0, if) = 0) 
D = 2m, F = 1.19m, A\ = A 2 = A 3 = Aa = 0. 
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Figure 17. Normalized \Ez\. Oblique incidence (0 — 5° , <p — 0, ip — 0) 
D = 2m, F = 1.19m, Ai = A 2 = A 3 = A 4 - 0. 
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Figure 18. 


Normalized |i£:n|. Oblique incidence (6 — 5° , 4> — — 0) 

D — 2m, F = 1.19m, A x = 0.1m, A 2 = -3, A z = A 4 = 0. 
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Normalized \E z\. Oblique incidence (6 = 5° , <j> = xp = 0). 
D = 2m, F — 1.19m, Ai — 0.1m, Ai — -3, A$ = A 4 = 0. 


Figure 20. 
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FORTRAN SUBROUTINE JOORAN 


SUBROUTINE JOORAN CALCULATES THE THREE COMPONENTS OF THE 
ELECTRIC FIELD SCATTERED FROM A PERFECT CONDUCTOR OF ARBI- 
TRARILY SHAPED SURFACE. THIS SURFACE IS SPECIFIED BY A 
SET OF TARGET POINTS. THE INCIDENT WAVE IS ASSUMED TO BE A 
PLANE WAVE WITH ORIENTATION DEFINED BY THREE EULERI AN 
ANGLES WITH RESPECT TO THE COORDINATE SYSTEM OF THE 
REFLECTOR. IT IS IMPORTANT TO MENTION THAT THE FIELD CAL- 
CULATED HERE IS NORMALIZED TO THE MAGNITUDE OF THE ELEC- 
TRIC FIELD OF THE INCIDENT WAVE. THE PLANE PERPENDICULAR 
TO THE DIRECTION OF PROPAGATION AND PASSES THROUGH THE 
ORIGIN IS THE ZERO PHASE PLANE. THE METHOD USED IS BASED 
ON A PHYSICAL OPTICS ANALYSIS. EDGE EFFECTS ARE INCLUDED 
IN THE SENSE THAT THE EFFECTS OF THE CHARGE AT THE BOUN- 
DARY OF THE ILLUMINATED SURFACE OF THE REFLECTOR HAVE BEEN 
COSIDERED (SEE SILVER "Microwave Antenna Theory and 
Design" section 5.8, Equations 62-69). FOR A DETAILED 
DESCRIPTION OF THE METHOD PLEASE, REFER TO THE REPORT 
ACCOMPANYING THIS SOURCE CODE. THIS CODE WAS DEVELOPED AT 
NORTH CAROLINA STATE UNIVERSITY BY DR. NICK E. BURIS AND 
DR. J. FRANK KAUFFMAN . SUPPORT FOR THIS STUDY WAS KINDLY 
OFFERED BY A NASA-LANGLEY RESEARCH CENTER GRANT. THIS CODE 
IS WRITTEN SO THAT IT CAN BE MODIFIED WITH THE MINIMUM 
POSSIBLE EFFORT. QUESTIONS ON THE OPERATION OF THIS CODE 
MAY BE ADDRESSED TO : 

Nick E. Buris 

University of Massachusetts 

Dept, of Electrical and Computer Eng. 

Marcus Hall Rm 20 

Amherst, MA 01003 

tel # (413) 545-2170 

ATTENTION : ALL LENGTHS IN METERS, FREQUENCY in MHz. 

SUBROUTINE JOORAN USES SUBROUTINE I B I RAN , 
THEREFORE, WHEN COMPILING MAKE SURE THAT 
YOU LINK TO A COMPILED VERSION OF IBIRAN. 

ARGUMENT DESCRIPTION 
INPUTS 

nd number of target points specifying the reflector 

( nd must be grater than 3 ) 


xdtil 

real array of dimension 
target points) 

nd ; 

(x 

coordinates 

of 

ydtil 

real array of dimension 
target points) 

nd ; 

(y 

coordinates 

of 

zdtil 

real array of dimension 
target points) 

nd ; 

( Z 

coordinates 

of 


nob number of observation points 

xotil real array of dimension nob; (x coordinates of 


A2 


yotil 

zotil 

aid 
a2d 
a3d 
f req 
ngauss 


Estil 


wk 

iwk 

xd 

yd 

zd 

xob 

yob 

zob 

x 

y 

z 

zx 

Ro 

SCAl 

SCA2 

SCA3 

gw 

gz 

u 

V 


observation points) 

real array of dimension nob; (y coordinates of 
observation points) 

real array of dimension nob; (z coordinates of 
observation points) 

Eulerian angle theta in degrees 

Eulerian angle phi in degrees 

Eulerian angle psi in degrees 

real variable ; frequency in MHz 

integer variable that determines the density of 

points used in the Gaussian Quadrature method of 

numerical integration. The larger ngauss the better 

the approximation. Available ngauss=4 , 6 , 8 , 10 , 12 , 16 , 48 . 

In the report, N and ng are used for ngauss. 

OUTPUT 

Complex array of dimension (3, nob). It contains the 
values of the scattered Electric field. 

Estil(j,k) = j-th component of E at k-th observ. point 
j = 1 => Ex, j = 2 => Ey, j = 3 => Ez 

WORK ARRAYS (OUTPUTS) 

real array of dimension 14*nd 

integer array of dimension 31*nd+l 

real array of dimension nd 

real array of dimension nd 

real array of dimension nd 

real array of dimension nob 

real array of dimension nob 

real array of dimension nob 

real array of dimension ( ngauss , ngauss ) 

real array of dimension ( ngauss , ngauss ) 

real array of dimension ( ngauss , ngauss ) 

real array of dimension ( ngauss , ngauss ) 

real array of dimension ( ngauss , ngauss ) 

complex array of dimension ( ngauss , ngauss ) 

complex array of dimension ( ngauss , ngauss ) 

complex array of dimension ( ngauss , ngauss ) 

double precision real array of dimension ngauss 

double precision real array of dimension ngauss 

real array of dimension ngauss 

real array of dimension ( ngauss , ngauss ) 


TO USE JOORAN SIMPLY CALL IT WITH A CALL STATEMENT LIKE 

call Jooran(nd,xdtil,ydtil,zdtil, nob ,xotil, yotil, zotil-, 

1 aid , a2d , a3d , f req , ngauss , Esti 1 , 

2 wk , iwk , xd , yd , zd , xob , yob , zob , x , y , z , zx , Ro , SCAl , SCA2 , SCA3 , 

3 gw,gz,u,v) 
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The first two lines of arguments contain the input argu- 
ments and the output argument Estil. The last two lines 
of arguments contain work arrays. The general user need 
not worry about these arrays. However, if changes in the 
program are needed for extensions and/or modifications, 
then the user should consult the theoretical part of this 
report and keep in mind that : 

1) xd,yd and zd are the coordinates of the target points in the 
ray system ( OXYZ ) . 

2) xob,yob and zob are the coordinates of the observation points 
in the ray system ( OXYZ ) . 

3) x,y,z,zx contain the values of x,y,z and Zx at the points defined 
by the grid ngauss X ngauss on each individual triangle. These 
values are used in the numerical integration scheme. 

4) Ro contains the values of the distance of the current observation 
point from each of the point of the grid mentioned in 3. 

5) SCAl,SCA2 and SCA3 contain some common factors that appear in 
the integral expressions for all the three components of the 
Electric field. Again, these values are calculated at the points 
of the grid mentioned in 3. 

6) gw,gz are the weights and the zeros of Legendre polynomials 
used in the Gaussian Quadrature integration. 

7) u,v are the coordinates of the points of the grid mentioned in 3 
expressed in the u-v coordinate system. In this coordinate system 
all the individual triangles are mapped onto a basic triangle. 
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As an example of using Jooran consider the following data file 


NASA87.DAT 


5000.000 frequency in MHz 

4 one of (4,6,8,10,12,16,48) 

5.000000 0 . 0000000E+00 0 . OOOOOOOE+OO Euler angles 

22 # of target points 

O.OOOOOOOOOOOOOOOOE+OO O.OOOOOOOOOOOOOOOOE+OO 0. 0000000000000000 E+00 
0.3333333333333333 0. 0000000 0000000 00E+00 2.1367522151277409E-02 

-0.1666666666666667 0.2886751345948129 2 . 1367522151277410E-02 

-0.1666666666666667 -0.2886751345948129 2 . 1367522151277409E-02 

0.6666666666666667 O.OOOOOOOOOOOOOOOOE+OO 8.5470088605109637E-02 

0.3333333333333333 0.5773502691896258 8 . 5470088605109636E-02 

-0.3333333333333333 0.5773502691896258 8 . 5470088605109639E-02 

-0.6666666666666667 -2 . 9 379182 5193 58776E-17 8 . 54700886 051 096 37E-0 2 

-0.3333333333333333 -0.5773502691896258 8 . 5470088605109637E-02 

0.3333333333333333 -0.5773502691896258 8 . 5470088605109639E-02 

1.000000000000000 O.OOOOOOOOOOOOOOOOE+OO 0.1923076993614967 

0.8660254037844387 0.5000000000000000 0.1923076993614967 

0.5000000000000000 0.8660254037844386 0.1923076993614967 

-2. 2034386889519082 E- 17 1.000000000000000 0.1923076993614967 

-0.5000000000000000 0.8660254037844387 0.1923076993614967 

-0.8660254037844386 0.5000000000000000 0.1923076993614967 

-1.000000000000000 -4. 4068773779038163 E- 17 0.1923076993614967 

-0.8660254037844387 -0.5000000000000000 0.1923076993614967 

-0.5000000000000000 -0.8660254037844386 0.1923076993614967 

-7. 2674717409587322 E- 17 -1.000000000000000 0.1923076993614967 

0.5000000000000000 -0.8660254037844387 0.1923076993614967 

0.8660254037844386 -0.5000000000000001 0.1923076993614967 

2 (nob # of observation points) 

O.OOOOOOOOOOOOOOOOE+OO O.OOOOOOOOOOOOOOOOE+OO 1.299999952316284 
4. 7746483236551285 E- 02 O.OOOOOOOOOOOOOOOOE+OO 1.299999952316284 
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This data file, nasa87.dat, is input to the following program : 


real xdtil(200) ,ydtil(200) , zdtil (200) 
real xotil(50) , yot i 1 ( 50 ) , zoti 1 ( 50 ) 
real xd (200), yd (200) ,zd (200), wk (14*200) 
real Xob( 50 ) , Yob( 50 ) , Zob( 50 ) 
complex Estil(3,50) 
integer nd , iwk ( 3 1 *2 0 0 + 1 * 1 ) 

real x ( 48 , 48 ) , y ( 48 , 48 ) , z ( 48 , 48 ) 
real zx ( 48 , 48 ) 
real Ro ( 48 , 48 ) 

complex SCAl (48,48) ,SCA2(48,48) 
complex SCA3(48,48) 

dimension gw (48) ,gz( 48) ,u(48) , v ( 48 , 48 ) 
double precision gw,gz 

open(l,file='nasa87.dat' ) 

read( 1 , * ) f req 

read( 1 , * ) ngauss 

read( 1 , * ) aid, a2d, a 3d 

read ( 1 , * ) nd 

do 1 k=l,nd 

read(l,*)xdtil(k) ,ydtil(k) , zdtil ( k ) 

1 continue 
read ( 1 , * ) nob 
do 2 k=l,nob 

read(l,*)xotil(k) ,yotil(k) , zotil ( k ) 

2 continue 

call Joor an ( nd,xd til, yd til, zdtil, nob, xotil,yo til, zotil, 

1 aid, a2d, a3d, f req, ngauss, Estil , 

2 wk,iwk,xd,yd,zd,xob,yob,zob,x,y,z,zx,Ro, SCAl , SCA2 , SCA3 , 

3 gw , gz , u , v ) 


open( 2, file=' nasout.dat' ) 


write ( 2 , * ) ' 
do 3 k=l,nob 
write ( 2 , * ) ' 

$ ' , ' , zotil ( k ) , 
wr i te ( 2 , * ) ' 
wr i te ( 2 , * ) ' 
wr i te ( 2 , * ) ' 

3 continue 
stop 
end 


,nob,' # of observation points' 

' , xoti 1 ( k ) , ' , ' , yoti 1 ( k ) , 

Xo til, Yotil, Zotil' 

' , Esti 1 ( 1 , k ) , 'Extil' 

' , Esti 1 ( 2 , k ) , ' Eyti 1 ' 

‘ , Estil ( 3 , k ) , ' Ez ti 1 ' 
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The output is sent to the file "nasout.dat” 


NASOUT.DAT 


2 # of observation points 

0 . OOOOOOOE+OO , 0 . OOOOOOOE+OO , 1.300000 

( -1.501289, -0.84 54 167) Ext il 
(-7.9167667E-06 ,1 . 51689 14E-04 ) Eytil 
( 0. 24 58099, -0.583 4 02 0)Ez til 
4.77 46483E-02 , 0 . 0000000E+00 , 1.300000 

( 0.62746 02, 0.181 089 4) Ext il 
(-1 . 3166843 E- 03, -4. 8433035E-04) Eytil 
( -0. 1880 506, 0.51 3989 5 )Ez til 


Xotil , Yotil , Zotil 


Xotil , Yotil , Zotil 


fields 


BUMP .1 


3 .2 .2 



0.0 0.4 0.8 1.2 1.6 


X (m) 


Figure 11. Normalized E. Observation points on the Ox axis. 

/ = 500 MHz, D = ?m, F = 1.19m, pi = 20, 

00;, = F, 6 = <j> — ip 0. Ai = 0.1m, Ai — -3, A3 — 0.2m, A. 
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APPENDIX B 


SUBROUTINE JOORAN FORTRAN CODE 
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C********************************************************************** 

THE FOLLOWING SUBROUTINE CALCULATES THE THREE COMPONENTS OF THE * 
ELECTRIC FIELD SCATTERED FROM AN ARBITRARY PERFECT CONDUCTOR . * 

THE INCIDENT WAVE IS ASSUMED TO BE A PLANE WAVE WITH ORIENTATION* 
DEFINED BY THREE EULERIAN ANGLES WITH RESPECT TO THE OPTICAL * 
AXIS OF THE REFLECTOR. IT IS IMPORTANT TO NOTE THAT THE FIELDS * 
ARE NORMALIZED TO THE MAGNITUDE OF THE ELECTRIC FIELD OF THE * 
INCIDENT WAVE. THE PLANE PERPENDICULAR TO THE DIRECTION OF PRO-* 
PAGATION AND PASSES THROUGH THE ORIGIN IS THE ZERO PHASE PLANE. * 
THE METHOD USED IS BASED ON A PHYSICAL OPTICS ANALYSIS AND EDGE * 
EFFECTS ARE INCLUDED. FOR A DETAILED DESCRIPTION OF THE METHOD,* 
PLEASE REFER TO THE REPORT ACCOMPANYING THIS SOURCE CODE . * 

THIS CODE WAS DEVELOPED AT NORTH CAROLINA STATE UNIVERSITY BY * 
DR. NICK E. BURIS AND DR. J. FRANK KAUFFMAN . SUPPORT FOR * 
THIS STUDY WAS KINDLY OFFERED BY A NASA-LANGLEY RESEARCH CENTER * 
GRANT . * 

THIS CODE IS WRITTEN SO THAT IT CAN BE MODIFIED WITH THE MINIMUM* 
POSSIBLE EFFORT. * 

★ 
★ 

ARGUMENT DESCRIPTION * 

INPUTS * 

* 

nd number of target points defining the reflector * 

xdtil real array of dimension nd; (x coordinates of * 

target points) * 

ydtil real array of dimension nd; (y coordinates of * 

target points) * 

zdtil real array of dimension nd; (z coordinates of * 

target points) * 

nob number of observation points * 

xotil real array of dimension nob; (x coordinates of * 

observation points) * 

yotil real array of dimension nob; (y coordinates of * 

observation points) * 

zotil real array of dimension nob; (z coordinates of * 

observation points) * 

aid Eulerian angle theta in degrees * 

a2d Eulerian angle phi in degrees * 

a3d Eulerian angle alpha in degrees * 

freq real variable ; frequency in MHz * 

ngauss integer variable that determines the density of * 

points used in the Gaussian Quadrature method of * 

numerical integration. The larger ngauss the better * 

the approximation. Available ngauss=4 , 6 , 8 , 10 , 12 , 16 , 48 . * 

OUTPUT * 

Estil Complex array of dimension (3, nob). It contains the * 

values of the scattered Electric field. * 

Estil(j,k) = j-th component of E at k-th observ. point* 

j = 1 => Ex, j = 2 => Ey, j = 3 => Ez * 

WORK ARRAYS * 

wk real array of dimension 14*nd * 
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c 

iwk 

integer array 

of dimension 31*nd+l 

* 

c 

xd 

real array of 

dimension nd 

★ 

c 

yd 

real array of 

dimension nd 

★ 

c 

zd 

real array of 

dimension nd 

* 

c 

xob 

real array of 

dimension nob 

* 

c 

yob 

real array of 

dimension nob 

* 

c 

zob 

real array of 

dimension nob 

* 

c 

X 

real array of 

dimension ( ngauss , ngauss ) 

ie 

c 

Y 

real array of 

dimension { ngauss , ngauss ) 

* 

c 

z 

real array of 

dimension ( ngauss , ngauss ) 

* 

c 

zx 

real array of 

dimension ( ngauss , ngauss ) 

★ 

c 

Ro 

real array of 

dimension ( ngauss , ngauss ) 

★ 

c 

SCAl 

complex array 

of dimension ( ngauss , ngauss ) 

★ 

c 

SCA2 

complex array 

of dimension ( ngauss , ngauss ) 

* 

c 

SCA3 

complex array 

of dimension ( ngauss , ngauss ) 

★ 

c 

gw 

double precision real array of dimension ngauss 

★ 

c 

gz 

double precision real array of dimension ngauss 

* 

c 

u 

real array of 

dimension ngauss 

* 

c 

V 

real array of 

dimension ( ngauss , ngauss ) 

* 

* 

c 

c* 

subroutine Jooran( nd,xdtil ,ydtil , zdtil ,nob,xotil ,yotil , zotil , 

1 aid , a2d , a3d , f req , ngauss , Esti 1 , 

2 wk , iwk , xd ,yd, zd, xob, yob, zob , x , y , z , zx , Ro , SCAl , SCA2 , SCA3 , 


3 gw,gz,u,v) 

real xdtil(nd) ,ydtil(nd) , zdtil ( nd ) 

real xotil(nob) , yot il ( nob ) ,zotil(nob) 

real xd(nd) , yd(nd) , zd(nd) ,xi ,yi , zi ,wk( 14*nd) 

real Xob(nob) , Yob(nob) ,Zob(nob) 

complex Estil(3,nob) 

integer nd, iop( 2 ) , iwk( 31*nd+l*l ) , ier r , nxi , nyi 


real x ( ngauss , ngauss ) ,y ( ngauss , ngauss ) , z ( ngauss , ngauss ) 
real zx ( ngauss , ngauss ) 
real Ro ( ngauss , ngauss ) 

complex SCAl ( ngauss , ngauss ) , SCA2 ( ngauss , ngauss ) 

complex SCA3 (ngauss, ngauss) 

dimension gw(ngauss) ,gz(ngauss) ,u(ngauss) ,v(ngauss, ngauss) 


double 

common 

COMMON 

1 


precision gw,gz 
/observ/ Xo , Yo , Zo , beta 

/IBCDPT/ X0,Y0,AP,BP,CP,DP,P00,Pl0,P20,P30,P40,P50, 
P01,PU,P21,P31,P41,P02,P12,P22,P32,P03,P1 
P23,P04,P14,P05, ITPV 


3 


/ 


common /vertex/ xl , x2 , x3 , yl , y2 , y3 
complex JKRZ 

complex Ex, Ey , Ez , Ex tot , Eytot , Eztot , Exknt , Eyknt , Ezknt 
external Ex,Ey,Ez 


pi=acos ( -1 . 0 ) 

Convert the angles from degrees into radians 
al r=ald*-pi/180 . 
a2r=a2d*pi/180 . 


C 


onoon o nonnoo oo no o ooo 
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a3 r=a3d*pi/180 . 

C Now calculate the elements of the Rotation Transformation Matrix, Rij 
cl=cos ( al r ) 
c2 = cos ( a2 r ) 
c3 = cos ( a3 r ) 
sl=s in ( al r ) 
s2=sin ( a2 r ) 
s3=sin( a3r ) 

Rll=c3*c2*cl-s3*s2 

Rl2=c3*s2*cl+s3*c2 

Rl3=-c3*sl 

R21=-s3*c2*cl-c3*s2 

R22=-s3*s2*cl+c3*c2 

R23=s3*sl 

R31=sl *c3 

R32=sl*s2 

R33=cl 

r r t t t t t t f t r t t r t t t r t t r t t t t t t r t t t t r r t r t t t t t r t t r r r t t r t r r r t r t t r r t r t t r t r t j 

Now Enter the Wavenumber ] 

make sure that f is in MHz ] 

beta=(2.*pi*freq)/300. 


iiii 'Now Enter the data points 


I t t t t t t r t i f t t t f r t t r t t f r t r t f 


] 

] 


call 


rotate(xdtil,ydtil,zdtil,xd,yd,zd,nd / alr,a2r,a3r) 


] 

] 




Now Enter The Observation Points 

read( 5 , * ) nob 
do 3 k=l,nob 

read( 5 , * ) Xob( k ) , Yob( k ) , Zob( k ) 

3 continue 

call rotate(xotil,yotil,zotil,xob, yob , zob , nob ,alr,a2r,a3r) 


] 

] 

1 

] 

] 

] 


t t t t t r j 

] 
] 
] 

(- 1 , 1 ) . 

i f ( ngauss . eq . 4 ) call gauss4 ( ngauss , gz , gw ) 
if ( ngauss . eq . 6 ) call gauss6 ( ngauss , gz , gw) 
if (ngauss. eq. 8) call gauss8 ( ngauss , gz , gw) 


r r r i t f i t r r > i t i i i i > t r i r t i t i i t r r r r r r t r r r r r r r i i i r r t / i > r t i r r t r r t r i 


create the (r,s) pairs. These are the points where the 
value of the integrand is needed in the r-s system. In 
this system the triangular domain is the square (-1,1) X 


noo nooo ooo no nonon 
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C 


if ( ngauss . eq . 10 ) 
i f ( ngauss . eq . 12 ) 
i f ( ngauss . eq . 16 ) 
i f ( ngauss . eq . 48 ) 


call gauslO ( ngauss , gz , gw) 
call gausl2 ( ngauss , gz , gw) 
call gausl6 ( ngauss , gz , gw ) 
call gaus48 ( ngauss , gz , gw) 


rtrftftttfttrrttftttrttttrrtfttttttrttirrrrtifrttittrttfftttftttttttt j 

Create the [u,v] pairs. These are the points where the ] 
value of the integrand is needed in the u-v system. In this ] 
system the domain is a right triangular defined by the points : ] 
( xl , yl ) <-> [0,0] (x2,y2) <-> [1,0] (x3,y3) <-> [0,1] ] 


call 


rstouv( ngauss , gz , u , v 


] 

] 


Now run ibiran for the first time 

Create the point mess 

iop( 1 ) =0 

iop ( 2 ) =0 

maxnxi=l 

nxi = l 

nyi = l 

xi- ( xd ( 1 ) +xd ( 2 ) +xd ( 3 ) +xd( 4 ) )/4 . 
yi = ( yd( 1 )+yd( 2 )+yd( 3 )+yd('4 ) )/4 . 


] 

] 

] 


call ibiran(nd,xd,yd,zd,iop, maxnxi ,nxi,xi,nyi,yi,zi,iwk,wk,ierr) 


nt=iwk ( 5 ) 
nl=iwk ( 6 ) 

write(2,*)"# of triangles =",nt," nl =",nl 


] 

] 


f / r t 9 t t 9 9 t f t t t r t f r t f t f r t r t r r r r t t t r t r t t r r t r t r t t t f r t r r t t r r r r t t 


This is the loop for each observation point ] 

do 6 lobe=l,nob 
Xo=Xob( lobe ) 

Yo=Yob( lobe ) 

Zo=Zob( lobe ) 


C 

C 

C 


Extot=( 0 .0,0.0) 
Eytot= ( 0 . 0 , 0 . 0 ) 
Eztot= ( 0 . 0 , 0 . 0 ) 


t t t t t r r t t t t t t t t r t t f t t r r t t f r t t t t r r t t t t t r r t r r t t t r t r t r t r f r 

This is the LOOP that deals with each triangle 


non oonon no noon 
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The barycentric point is found of each triangle and ] 
the ibiran is called in order to calculate the coef- ] 
ficients of the polynomial in the u-v system at the ] 
triangle in question. ] 

do 5 knt=l,nt 

iop( 1 ) =0 

iop( 2 )=2 

nxi=*l 

nyi = l 

maxnxi=l 

kk=15+3*knt 


xl=>xd( iwk ( kk-2 ) ) 
yl=yd ( iwk ( kk-2 ) ) 
x2»xd ( iwk ( kk-1 ) ) 
y2-yd ( iwk ( kk-1 ) ) 
x3*xd ( iwk ( kk ) ) 
y3*yd( iwk ( kk ) ) 


a=x2-x 1 
b*x3-xl 
c=*y2-yl 
d=y 3-y 1 

abcd= a*d - b*c 

UparX= d / abed 
VparX« -c / abed 

The coordinates of the Barycentric point of the knt-th ] 
triangle follow ] 

xbaryc = (xl + x2 + x3) / 3.0 
ybaryc = (yl + y2 + y3 ) / 3.0 


Run IBIRAN in order to create the coefficients of the poly- ] 
nomial that describes the surface in the u-v system. To make] 
sure that the correct triangle is considered IBIRAN is asked ] 
to run for the BARYCENTRIC POINT of the triangle. This point] 
lies in the triangle no matter what !!!!!!!!!!!!!!!!!!!!!!! ] 

call ibiran( nd , xd,yd , zd, iop,maxnxi , nxi , xbaryc , nyi , ybaryc , 

1 zbaryc , iwk , wk , ier r ) 


] 

Create the ( x , y ) paTFsT These are the points where the ] 

value of the integrand is needed in the x-y system. In ] 


oonoonn onoooo non 
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this system the triangular domain is defined by the points 
( xl , yl ) , (x2,y2) and (x3,y3) 

call uvtoxy(ngauss,x,y,u,v,xl,x2,x3,yl,y2,y3) 

Create the Z and Zx values. Z(k,j) is the value of 
the reflector surface Z=Z(x,y) calculated at the point 
x=x ( k , j ) , y=y(k,j) . Zx(k,j) is the value of the partial 

derivative of Z with respect to x at the aforementioned 
point - ( x ( k , j ) , y(k,j) 


] 

] 

] 

] 

] 

] 

] 

] 

] 


49 

50 


do 50 k=l,ngauss 
uk=u ( k ) 

do 49 j=l,ngauss 
vk j=v ( k , j ) 

as0 = P00+ vk j * ( P01+ vk j * ( P02+ vkj*(P03+ vkj*(P04+ vkj*P05)))) 

asl=PlO+ vk j * ( P11+ vk j * ( P12+ vkj*(Pl3+ vkj*Pl4))) 

as2 = P20+ vk j * ( P21+ vkj*(P22+ vkj*P23)) 

as3 = P30+ vk j * ( P31+ vkj*P32) 

as4=P40+ vkj*P41 

Z(k,j)=as0+ uk*(asl+ uk*(as2+ uk*(as3+ uk*(as4+ uk*P50)))) 
Zu=asl+ uk*(2.*as2+ uk*(3.*as3+ uk*(4.*as4+ uk*5.*P50))) 
as0=P01+ vkj*(2.*P02+ vkj*(3.*P03+ vkj*(4.*P04+ vk j *5 . *P05 ) ) ) 
asl = Pll+ vk j * ( 2 . *P12+ vk j * ( 3 . *Pl 3 + vkj*4.*Pl4)) 
as2 = P21+ vk j * ( 2 . *P22+ vkj*3.*P23) 
as3=P31+ vkj*2.*P32 

Zv=as0+ uk*(asl+ uk*(as2+ uk*(as3+ uk*P41 ))) 

Zx(k,j)= Zu*UparX + Zv*Vparx 

continue 

continue 


Create the array Ro(k,j). This array contains the distance 
of the point ( x ( k , j ) ,y ( k , j ) , z ( k , j ) ) from the observation 

point ( Xo , Yo , Z o ) 

ALSO 

Create the complex arrays SCA(k,j). These arrays contain ] 
the common ceofficients in front of all the field compo- ] 
nents. ] 

do 58 k=l,ngauss 
do 57 j=l,ngauss 

Ro ( k , j ) =sqr t ( ( x ( k , j ) -Xo ) * ( x ( k , j ) -Xo ) + 

1 ( y ( k , j )-Yo ) * ( y ( k , j ) -Yo ) + 

2 ( z( k, j )-Zo ) * ( z( k, j )-Zo) ) 

RR=Ro ( k , j ) 

JKRZ=cmplx ( 0.0, -beta * (RR-z(k,j)) ) 

SCAl ( k , j ) = cexp(JKRZ) / RR 

SCA2 ( k , j ) = cmplx(beta*beta - 3.0/(RR*RR), -3 . 0 *beta/RR ) * 

1 ( Xo-x ( k , j ) + ( Zo-z ( k , j ) ) *Zx( k , j ) )/RR 

SCA3 ( k , j ) = cmplx ( 1 . 0/( RR*RR ) -beta*beta , beta/RR ) 


no n onoonnnooo no no no 
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57 continue 

58 continue 


Call the "CALINT" subroutine to calculate the 
Integrals on each triangular domain 


call calint(Ex,ngauss,x,y,z,zx,Ro, SCAl , SCA2 ,SCA3,Exknt,gw,gz,u,v) 
call cal int ( Ey , ngauss , x , y , z, zx,Ro,SCAl , SCA2 , SCA3 , Eyknt , gw , gz , u , v ) 
call cal int ( Ez , ngauss , x , y , z , zx , Ro , SCAl , SCA2 , SCA3 , Ezknt , gw, gz , u , v ) 


* 


Extot=Extot + Exknt 
Eytot=Eytot + Eyknt 
Eztot=Eztot + Ezknt 
5 continue 

Extot=Extot/cmplx ( 1 . 0 , 2.0*pi*beta) 
Eytot=Eytot/cmplx ( 1 . 0 , 2.0*pi*beta) 
Eztot=Eztot/cmplx ( 1 . 0 , 2.0*pi*beta) 


" " Transform the E-fields back to the tilded coordinate system '''''] 

] 


Estil ( 1 , lobe )=Rll *Extot+R21 *Eytot+R3 l*Ez tot 
Estil ( 2, lobe )=Rl2*Extot+R22*Eytot+R32*Eztot 
Estil ( 3 , lobe )=Rl3*Extot+R23*Eytot+R33*Eztot 


write( 3 , * ) ' 
write ( 3 , * ) ' 
write ( 3 , * ) ' 
write ( 3 , * ) ' 
wr ite ( 3 , * ) 1 
wri te ( 3 , * ) ' 
wr i te ( 3 , * ) ' 
wr i te ( 3 , * ) ' 


" ,Rll*Xo+R21*Yo+R31*Zo, "correct" 

" , Rl2*Xo+R22 *Yo+R32 *Zo , "correct" 

" , Rl 3 *Xo+R2 3 * Yo+R3 3 *Zo , " correct" 

" ,Xotil( lobe) , Yoti 1 ( lobe ) , Zotil(lobe) , "coord. 
", Estil ( 1 , lobe ), " Esxtil" 

", Estil ( 2 , lobe ), " Esytil" 

", Estil ( 3 , lobe ), " Esztil" 


72 

in til" 


] 


6 continue 

return 

end 


72 


This subroutine rotates the antenna as well as the observation 
points into the incident wave coordinate system ! 


I 
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subroutine rotate ( Axt , Ayt , Azt , Ax , Ay , Az , n , al , a 2 , a 3 ) 
dimension Ax ( n ) , Ay ( n ) , Az ( n ) , Axt ( n ) , Ayt ( n ) , Azt ( n ) 

cl=cos(al) 
c2=cos ( a2 ) 
c3=cos ( a3 ) 
sl=sin(al) 
s2=sin ( a2 ) 
s3=sin( a3 ) 

Rll=c3*c2*cl-s3*s2 

Rl2=c3*s2*cl+s3*c2 

Rl3=-c3*sl 

R21=-s3*c2*cl-c3*s2 

R22=-s3*s2*cl+c3*c2 

R23=s3*sl 

R31=sl*c3 

R32=sl*s2 

R33=cl 

do 1 i=l,n 

Ax ( i ) = Rll*Axt(i) + Rl2*Ayt ( i ) + Rl3*Azt(i) 

Ay ( i ) =» R21*Axt(i) + R22*Ayt ( i ) + R23*Azt(i) 

Az(i) = R31*Axt(i) + R32*Ayt(i) + R33*Azt(i) 

1 continue 

return 

end 


C 72 

subroutine calint(f,n,x,y,z,zx,Ro, SCAl , SCA2 , SCA3, result, gw, gz,u,v) 
dimension u(n) ,v(n,n) ,x(n,n) ,y(n,n) ,z(n,n) ,zx(n,n) 
real Ro(n,n) 

complex SCAl ( n , n ) , SCA2 ( n , n ) , SCA3 ( n , n ) 
dimension gw(n),gz(n) 
double precision gw,gz 

COMMON /IBCDPT/ XO , YO , AP , BP , CP , DP , POO , P10 , P20 , P30 , P40 , P50 , 

1 P01 / P11,P21,P31,P41,P02,P12,P22,P32,P03,P13, 

2 P23,P04,P14,P05, ITPV 
common /vertex/ xl , x2 , x3 , yl , y2 , y3 
complex f , Tkj , result , res 

external f 


This program calculates the double integral 

I = f(x,y,z,zx) dxdy on a triangular domain 

where z=z(x,y) and zx is partialz/partialx (also function of x,y). 
The domain is determined by the vertices of the triangle 
( xl , yl ) , (x2,y2) and (x3,y3) GIVEN COUNTERCLOCKWISE 1!! 

If the vertices are not arranged counterclockwise then 


i 



oooo no nooooonnooo 
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the value of I is the opposite of what it should be It!! 

abed as calculated below is the Jacobian of the mapping 
from the x-y to the u-v system. It is equal to the ratio 
of the areas of the triangular domain in the two systems. 
Since in the u-v system the triangle is always given by 
the points [0,0], [1,0] and [0,1], in this system its area 
is 1/2 . Thus, it is not a surprise that abed equals twice 
the area of the triangle in the x-y system (defined by the 
points (xl,yl), (x2,y2) and (x3,y3) ) !!!!!!!!!!!!!!!!!!!!! 


a=x2-xl 
b=x3-xl 
c=y2-yl 
d=y3-yl 

abcd= a*d - b*c 
i f ( abed . eq . 0 . 0 ) goto 666 

Calculate the integral by using Gaussian Quadrature 
in 2-dimensions. For the algorithm see my notes. 

This part of the subroutine calculates the values of the function 
T(r,s) at the points that are used by the Gauss-Legendre 
quadrature method. These values are used in the main 
in a summation scheme to provide the approximation to the 
c integral. Since in this program we have used m=n, r(k) 

c and s ( k ) are exactly the same. For a more general approach 

c read my notes, 

i c 

I c A very important note should be made here. abed is the 

c value of the Jacobian of the mapping from the x-y to u-v. 

c Since abed is not the absolute value of the Jacobian one needs 

c to make sure that the points (xl,yl), (x2,y2) & (x3,y3) are 

c ordered counterclockwise. For more details look up my notes. 


result=( 0 . 0 , 0 . 0 ) 
res=(0.0,0.0) 
do 2 k=l , n 
do 1 j=l,n 

Tk j=0 . 125 * (l.-gz(k)) * abed * 

1 f( x( k, j ) ,y(k, j ) ,z( k, j ) , zx(k, j ) ,Ro(k, j ) , 

2 SCAl ( k , j ) , SCA2 ( k , j ) , SCA3 ( k , j ) ) 

res=res + gw(j) * Tkj 

1 continue 

result=result + gw(k) * res 
res= ( 0 . 0 , 0 . 0 ) 

2 continue 
go to 667 

666 write ( 0 ,*) "There is something wrong with the triangle.", 
l"Colinear vertices." 

667 continue 



n o 


return 

end 


BIO 


subroutine uvtoxy ( n,x,y,u,v,xl,x2,x3,yl,y2,y3) 
dimension u(n) ,v(n,n) ,x(n,n) ,y(n,n) 


a=x2-xl 
b=*x3-xl 
c=y2-yl 
d=y 3-yl 

abcd= a*d - b*c 

do 2 k = l , n 
do 1 j-l,n 

x(k,j)=xl + a * u ( k ) + b * v(k,j) 
y(k,j)=yl + c * u(k) + d * v(k,j) 

1 continue 

2 continue 

return 

end 


subroutine r stouv ( n , gz , u , v ) 
dimension gz ( n ) , u ( n ) , v ( n , n ) 
double precision gz 

do 1 k=l , n 

u( k ) =0 . 5 + 0 . 5*gz ( k ) 

1 continue 

do 3 k=l , n 
do 2 j=l,n 

v ( k , j ) = ( 1 . 0-gz ( k ) ) * ( 1 . 0+gz ( j ) ) / 4.0 

2 continue 

3 continue 

return 

end 


complex function Ex ( x , y , z , zx , Ro , SCAl , SCA2 , SCA3 ) 

Here Ro and the SCA's are constants since Ex is called by ] 
x(k,j),y(k,j),z(k,j),zx(k,j),Ro(k,j),SCA's(k,j) ] 

common /observ/ Xo , Yo , Zo , beta 
complex SCAl , SCA2 , SCA3 


I 


o o no 


Ex=SCAl * ( SCA2* (Xo-x )/Ro + SCA3 ) 

return 

end 
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complex function Ey ( x , y , z , zx , Ro , SCAl , SCA2 , SCA3 ) 

Here Ro and the SCA's are constants since Ey is called by ] 
x(k,j),y(k,j),z(k,j),zx(k,j),Ro(k,j),SCA's(k,j) j 

common /observ/ Xo , Yo , Zo , beta 
complex SCAl , SCA2 , SCA3 

Ey=SCAl * SCA2* ( Yo-y )/Ro 

return 

end 


complex function Ez ( x , y , z , zx , Ro , SCAl , SCA2 , SCA3 ) 

Here Ro and the SCA's are constants since Ez is called by ] 
x(k,j ) ,y(k,j ) ,z(k,j ) , zx ( k , j ) ,Ro(k,j ) , SCA' s ( k , j ) ] 

common /observ/ Xo , Yo , Zo , beta 
complex SCAl , SCA2 , SCA3 

Ez=SCAl * ( SCA2* ( Zo-z )/Ro + SCA3*zx ) 

return 

end 


subroutine gauss4 ( n, z ,w) 
double precision z,w 
dimension z(n),w(n) 
n=4 

nn=n/2 

z(l)=*.339981043584856d0 
z( 2)=.861136311594053d0 

w(l)=.652145154862546d0 

w(2)=.347854845137454d0 

do 1 k®l,nn 
z ( nn+k ) =-z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l,nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 

2 continue 
return 
end 


c 
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c 

subroutine gauss6 ( n , z , w ) 
double precision z,w 
dimension z(n),w(n) 
n=6 

nn=n/2 

z(l)=.238619186083197d0 
z(2)=.661209386466265d0 
z( 3)=.932469514203152d0 

w ( 1 )«. 467913934572691d0 
w ( 2 )=. 360761573048139d0 
w ( 3)=.171324492379170d0 

do 1 k=l,nn 
z ( nn+k ) =-z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l,nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 

2 continue 
return 
end 

c 

c 

subroutine gauss8 ( n, z , w) 
double precision z,w 
dimension z(8),w(8) 


c 

c The array z contains the zeros of the 

c of the 8th degree. The array w conta 

c weights for the Gaussian quadrature, 

c are given by 

c 

c w( i ) = 2 / ( ( 1-z ( i ) * * 2 ) * ( P' 

c 

c where P(z(i))=0 i=l,2,...,8 

c and P' ( z )= dP( z )/dz 

c 


n=8 

nn=n/2 

z(l)=.183434642495650d0 
z( 2 )=. 525532409916329d0 
z(3)=.796666477413627d0 
z(4)=.960289856497536d0 


Legendre polynomial 
ns the appropriate 
The values of w(k) 


z ( i ) ) **2 ) ) 


o o 
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w(l)=.362683783378362d0 
w(2)=.313706645877887d0 
w( 3)=.222381034453374d0 
w( 4 )=. 10122853629037 6d0 


do 1 k=l,nn 
z ( nn+k ) =-z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l,nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 

2 continue 
return 
end 


subroutine gausl 0 ( n , z , w ) 
double precision z,w 
dimension z(n),w(n) 

n=10 

nn=n/2 

z(l)=.148874338981631d0 
z(2)=.433395394129247d0 
z(3)=.679409568299024d0 
z( 4)=.865063366688985d0 
z( 5)=.973906528517172d0 

w(l)=.295524224714753d0 
w ( 2 )=. 269266719 3 09996d0 
w ( 3 )=. 219086362515982d0 
w(4)=.149451349150581d0 
w(5)=.066671344308688d0 

do 1 k = 1 , nn 
z ( nn+k ) =»— z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l,nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 
continue 
return 
end 


2 


no on 


Bl 4 


subroutine gausl2 ( n , z , w ) 
double precision z,w 
dimension z(n),w(n) 
n-12 
nn=n/2 

z(l)=. 12523340 8 5 11469d0 
z(2)=.367831498998180d0 
z( 3 ) = . 587317954286617d0 
z(4)=.769902674194305d0 
z(5)=.904117256370475d0 
z( 6 )=.981560634246719d0 

w(l)=. 24914704581340 3d0 
w(2)=.233492536538355d0 
w( 3 ) = . 203167426723066d0 
w ( 4 ) = .160078328543346d0 
w( 5)=.106939325995318d0 
w(6)=.047175336386512d0 

do 1 k=l,nn 
z ( nn+k ) =-z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l , nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 

2 continue 
return 
end 


subroutine gausl 6 ( n , z , w ) 
double precision z,w 
dimension z(8),w(8) 


c 

c The array z contains the zeros of the Legendre polynomial 

c of the 16th degree. The array w contains the appropriate 

c weights for the Gaussian quadrature. The values of w(k) 

c are given by 

c 

c w(i)= 2 / ( ( 1-z ( i ) **2 ) * ( P'(z(i))**2 ) ) 

c 

c where P(z(i))=0 i=l,2,...,16 

c and P' ( z ) = dP(z)/dz 

c 


n=16 
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nn=n/2 

z(l)=. 09501250983663744018 5d0 
z(2)=.281603550779258913230d0 
z(3)=. 4 58 01677765722738634 2d0 
z( 4 )=. 61787624440264374844 7d0 
z(5)=. 755404408 35500303389 5d0 
z( 6 )=. 86 56 312023878317 4388 OdO 
z( 7 )=. 94457502307323257607 8d0 
z( 8)=.989400934991649932596d0 


w( 1 ) = . 1894 5061 0455068 496285d0 
w( 2 )=. 18260341 50 44923558867d0 
w ( 3)=.169156519395002538189d0 
w(4)=.149595988816576732081d0 
w( 5) =.12 46289712 5553 38720 52d0 
w( 6 )=. 09 51 58 511682 49278 481 OdO 
w(7)=. 06225352393864789286 3d0 
w( 8 )=.027152459411754094852d0 


do 1 k=l,nn 
z ( nn+k ) =-z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l,nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 

2 continue 


return 

end 


subroutine gaus48 ( n, z , w) 
c 

c The array z contains the zeros of the Legendre polynomial 

c of the 48th degree. The array w contains the appropriate 

c weights for the Gaussian quadrature. The values of w(k) 

c are given by 

c 

c w(i)= 2 / ( ( 1-z ( i ) **2 ) * ( P'(z(i))**2 ) ) 

c 

c where P(z(i))=0 i=l,2,...,48 

c and P'(z)= dP(z)/dz 

c 

dimension z(n),w(n) 
i double precision z,w 
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n=48 

nn=n/2 

z ( 1 ) = . 032380170962869362033d0 
z ( 2 ) =.097004699209462698930d0 
z ( 3 ) =. 161222356068891718056d0 
z ( 4 ) =.224763790394689061225d0 
z ( 5 ) =.287362487355455576736d0 
z ( 6 ) =. 348755886292160738160d0 
z ( 7 ) =.408686481990716729916d0 
z ( 8 ) =.466902904750958404545d0 
z ( 9 ) =. 523160974722233033678d0 
z(10)=. 577224726 08 3972703818d0 
z(ll)=.628687396776513623995d0 
z( 12) =.6778723796 3266 3905212d0 
z(13)=. 72403413092 38 1465467 4d0 
z(14)=.767159032515740339254d0 
z( 15)=.807066204029442627083d0 
z(16)=.843588261624393530711d0 
z( 17)=. 87657202027424788590 6d0 
z( 18) =.9 058791 3671 5569672822d0 
z( 19)=. 93138669070655433311 4d0 
z( 20) =.9529877031604 3086 072 3d0 
z( 21) =.97059159 2 5462 47250461d0 
z( 22) =.984124 58 37228268 5774 5d0 
z(23)=. 99 35301722663 507 57 548d0 
z( 24)=. 99877 10072 52 426118601d0 


w ( 1 ) =. 064737696812683922503d0 
w(2) =.064466164435950082207d0 
w( 3 ) =. 063924238584648186624d0 
w( 4 ) =. 063114192286254025657d0 
w( 5 ) =.062039423159892663904d0 
w( 6 ) =.060704439165893880053d0 
w( 7 ) =. 059114839698395635746d0 
w( 8 ) =.057277292100403215705d0 
w( 9 ) =.055199503699984162868d0 
w ( 10)=. 05289018948 519 3667096d0 
w(ll)=.050359035553854474958d0 
w ( 12)=. 0476166 58 49249047 4 826d0 
w( 13)=. 044674 5608 56694280419d0 
w( 14)=. 04154508294346474921 4d0 
w( 15)=. 03824135106583070631 7d0 
w( 16 ) = . 0 3 477722256 477 0 4 3889 3d0 
w( 17)=. 03116722783279808890 2d0 
w( 18)=.027426509708356948200d0 
w( 19)=. 023 57 076083932437914 ldO 
w( 20)=. 0196161604 57 355527 81 4d0 
w( 21)=. 015 57 931 57 2294384 87 28d0 
w( 22) = . 011477234 57923453949 OdO 
w( 2 3 ) = . 007 3 27 553901276 2621 02d0 
w( 24)=. 0031533460 52 30 58 3863 3d0 
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do 1 k=l,nn 
z ( nn+k ) --z ( k ) 
w( nn+k ) =w( k ) 

1 continue 

do 2 k=l ; nn 
z ( k ) =-z ( n+l-k ) 
w( k ) =w( n+l-k ) 

2 continue 
return 
end 


c 


